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Abstract 

Inputs to signaling pathways can have complex statistics that depend on the environment and on the behavioral response 
to previous stimuli. Such behavioral feedback is particularly important in navigation. Successful navigation relies on proper 
coupling between sensors, which gather information during motion, and actuators, which control behavior. Because 
reorientation conditions future inputs, behavioral feedback can place sensors and actuators in an operational regime 
different from the resting state. How then can organisms maintain proper information transfer through the pathway while 
navigating diverse environments? In bacterial chemotaxis, robust performance is often attributed to the zero integral 
feedback control of the sensor, which guarantees that activity returns to resting state when the input remains constant. 
While this property provides sensitivity over a wide range of signal intensities, it remains unclear how other parameters such 
as adaptation rate and adapted activity affect chemotactic performance, especially when considering that the swimming 
behavior of the cell determines the input signal. We examine this issue using analytical models and simulations that 
incorporate recent experimental evidences about behavioral feedback and flagellar motor adaptation. By focusing on how 
sensory information carried by the response regulator is best utilized by the motor, we identify an operational regime that 
maximizes drift velocity along chemical concentration gradients for a wide range of environments and sensor adaptation 
rates. This optimal regime is outside the dynamic range of the motor response, but maximizes the contrast between run 
duration up and down gradients. In steep gradients, the feedback from chemotactic drift can push the system through a 
bifurcation. This creates a non-chemotactic state that traps cells unless the motor is allowed to adapt. Although motor 
adaptation helps, we find that as the strength of the feedback increases individual phenotypes cannot maintain the optimal 
operational regime in all environments, suggesting that diversity could be beneficial. 
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Introduction 

Escherichia coli cells navigate their environment by alternating 
straight runs with direction-changing tumbles to perform a 
random walk. During a run, the flagellar motors spin counter- 
clockwise (CCW) and propel the cell at constant speed in one 
direction, which changes slowly due to rotational diffusion. Runs 
are terminated when one or more motors start rotating clockwise 
(CW), which causes the cell to tumble [1—3]. Cells are able to bias 
their random walk toward favorable conditions using a two- 
component signal transduction pathway that detects changes in 
signal intensity during runs and modulates the probability to 
tumble accordingly, resulting in extended runs in the desired 
direction and net drift velocity in the direction of the gradient [4] . 

The sensory module of the chemotaxis pathway (Figure 1A) 
consists of large clusters of receptor proteins that bind signal 
molecules to modulate rapidly (<0. 1 s) the activity of an associated 
histidine kinase, CheA [5-7]. The high gain of the receptor cluster 
is coupled to negative integral feedback control [8-10], mediated 



by slow (~ 1-30 seconds) methylation and demethylation of the 
receptors by CheR and CheB, respectively [11—13]. This allows 
the receptors to adapt to a constant background signal while 
maintaining sensitivity over a wide range of concentrations 
[14,15]. For example, when cells are stimulated with a step of 
aspartate, the activity of the receptors returns nearly precisely to its 
pre-stimulus level after a transient response (Figure IB first line). 
While precise adaptation does not hold when receptors become 
saturated, adaptation with a precision above 80% has been 
measured for many relevant signals within the micromolar range 
[16]. Precise adaptation is an important feature of bacterial 
chemotaxis because it provides robustness by implementing a 
"maximin" strategy that guarantees at least minimum chemotactic 
performance in any environmental condition [17]. 

The activity of the sensory module is relayed through a 
diffusible response regulator CheY to the flagellar motors, which 
act as the actuator, (Figure 1A). When phosphorylated by CheA, 
CheY-P binds to the motor subunit FliM and increases the 
probability of the motor to switch from CCW to CW [18]. Fast 
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Author Summary 

The biased random walk is a fundamental strategy used by 
many organisms to navigate their environment. Drift along 
the desired direction is achieved by reducing the proba- 
bility to reorient whenever conditions improve. In the 
chemotaxis system of Escherichia coli, this is accomplished 
with a sensory module that implements negative integral 
feedback control, the output of which is relayed to the 
flagellar motors (the actuators) by a response regulator to 
control the probability to change direction. The proper 
dynamical coupling between sensor and actuator is critical 
for the performance of the random walker. Here, we 
identify an optimal regime for this coupling that maximiz- 
es drift velocity in the direction of the gradient in multiple 
environments. Our analysis reveals that feedback of the 
behavior onto the system in steep gradients can constrain 
individual cell performance, by causing bi-stable behavior 
that can trap cells in non-chemotactic states. These 
limitations are inherent in the biased random walk strategy 
with integral feedback control, but can be alleviated if the 
output of the pathway adapts, as recently characterized for 
the flagellar motors in Escherichia coli. 



dephosphorylation of CheY-P by the phosphatase CheZ ensures 
rapid transfer of information from the sensor to the actuators. 

The CW bias of the flagellar motor, which defines the tumbling 
probability [2], is a sensitive function of the CheY-P concentration 
(Hill coefficient >10, Figure 1A) [19,20]. The capability of the 
system to maintain the CheY-P concentration within the tight 
dynamic range of the motor CW bias response function (Figure 1 A) 
is often used to investigate robustness to fluctuations in protein 
concentrations and receptor activity [21,22]. An important 
underlying assumption is that performance is maximized when 
the motor converts small variations in CheY-P into large changes 
in CW bias. 

However, recent experiments and theory suggest that the 
coupling between sensor and actuator is more complex than 
previously thought. First, the flagellar motors partially adapt to 
persistent stimulus [23,24]. Second, the motor CW bias response 
to CheY-P is steeper than previously reported, further restricting 
the dynamic range of the motor response to CheY-P fluctuations 
[20] . Finally, in exponential ramps of chemoattractant, the CheY- 
P concentration reaches a dynamical equilibrium, T m , hereafter 
called operational CheY-P, distinct from the adapted CheY-P 
concentration, Tg, that the cell maintains in constant uniform 
environments [25,26] (Figure IB second line). For each of these 
three findings, the characterization of the internal dynamics of the 
signaling pathway was performed on immobilized cells using 
experimentally controlled input signals. However, for cells 
swimming freely in chemical gradients, the input signal dynamics 
are determined by the chemotactic response of the cell, creating a 
feedback of the behavior onto the input signal [27] (Figure IB 
third line). Because of this behavioral feedback, it remains unclear 
how the multiple time scales of the system, from signal detection to 
motor response, ultimately determine chemotactic drift perfor- 
mance. 

Here, we use analytical models and stochastic simulations of 
individual cells to examine the consequences of these new 
observations for our understanding of the bacterial chemotaxis 
strategy. Clonal populations of chemotactic E. coli grown in 
homogeneous conditions exhibit significant cell-to-cell phenotypic 
variability, with adaptation times ranging from 1 to 30 seconds 
[28-31], and motor clockwise bias ranging from 0.1 to 0.4 [3,31]. 



Therefore, we consider how different combinations of adaptation 
times and motor clockwise biases, which define a cell's phenotype, 
affect individual cell chemotactic drift velocity in different 
environments. In a phenotypically diverse population, different 
phenotypes of that population may perform best in different 
environments. 

Focusing on how information transfer from sensor to actuator 
affects chemotactic performance, we analyze the dynamical 
relationship between the operational regime of CheY-P, T m , and 
the drift velocity, Vo, as a function of the phenotype and gradient 
steepness. We show that there is a unique operational regime of 
the sensor with respect to the motor that maximizes drift velocity 
in the direction of the gradient by maximizing the contrast 
between runs up and down the gradient, and not by maximizing 
the CW bias response. We characterize the performance trade-off 
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Figure 1. Dynamical coupling between the sensor and the 
actuator in the bacterial chemotaxis system. A. The bacterial 
chemotaxis system is composed of a sensor module (receptor-kinase 
complexes) and an actuator module (flagellar motors) coupled through 
the phosphorylated form of CheY. Both modules are ultra-sensitive and 
adapt to their respective input signals. Maintaining the output of the 
sensor within the right range relative to the actuator is critical for 
chemotaxis performance. B. Diagrams of the CheY-P concentration 
response to different signals. First line: when cells are immobilized onto 
a slide, a step stimulus of attractant (e.g. methylaspartate) causes a 
sudden decrease in CheY-P concentration followed by a slower 
adaptation. Because of the negative integral feedback architecture of 
the sensor module, CheY-P adapts back to its pre-stimulus level, the 
adapted CheY-P concentration, Y 0 . Second line: when immobilized cells 
are exposed to an exponential ramp in time of the same stimulus, the 
system, which is log sensing, experiences a constant "force" and adapts 
towards an operational CheY-P concentration, Y m , lower than the 
adapted level Y 0 . This deviation of CheY-P activity from Y 0 to Y m 
changes the coupling between sensor and actuator. Third line: when 
cells are swimming in a gradient of attractant, their biased random walk 
causes them to climb the gradient. The average drift velocity of the cell 
up a chemical gradient affects the average input signal experienced by 
the cell. This creates a feedback of the behavior onto the input signal, 
which in turn can significantly affect the operating concentration of 
CheY-P and thus the coupling between sensor and actuator. 
doi:10.1371/journal.pcbi.1003694.g001 
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Figure 2. Simulated and theoretical drift velocity l/ D in exponential gradient of aspartate L 0 e^ x . A. V D as a function of the adapted CheY-P 
concentration Y 0 , in a shallow gradient (L 0 = 200 uM and =5,000 urn) for cells with adaptation times x = 5 (blue), 10 (green), and 30 seconds (red). 
V D is the average velocity of 10,000 identical cells between f = 60 and 300 seconds (dots: stochastic simulations; lines: analytical solution from Eq. (3); 
grey: motor CW bias response curve. B. Expected trajectories of CheY-P concentration Y(F(t)) for cells running in one dimension up (green) or down 
(red) in a gradient (integration of Eqs. (2) and (5), see Text SI; t = 30 s, g _1 =5,000 um, V(F,) = 2.4 uM and 3 uM). Expected run, /.^ (dotted line), and 
tumble, Xjg (dashed line), durations as a function of Y 0 . Expected run duration along a given direction T m = (2D r +i R0 )~' (solid black line) is limited by 
rotational diffusion (D r = 0.062 rad 2 s _1 ). Grey: motor CW bias. C. Same as A (t = 10 s) but with the rotational diffusion constant D r = 0.031 (red), 0.062 
(green), and 0.124 (blue) rad 2 s _1 . Dotted lines: expected run duration in a given direction. D. Same as A (t = 10 s) but with the motor switching rate 
co = 2.6 (red), 1.3 (green), and 0.65 (blue) s _1 . Dotted lines: expected run duration in a given direction. 
doi:1 0.1 371 /journal.pcbi.1 003694.g002 



faced by individual cells with different combinations of phenotypic 
parameters (such as, adapted CheY-P concentration, T 0 , receptor 
adaptation time, T, and cell resting tumble bias, TB ( ]). 

Results 

Maximizing contrast in run durations rather than CW bias 
response maximizes chemotactic drift velocity 

Previous studies have examined how E. coli chemotactic drift 
velocity along a one dimensional gradient depends on the 
adaptation time [25,32], the shape of the response function of 
the sensory module [17,33,34] and also behavioral feedback [27]. 
Instead, we first focus on how the coupling between the sensors 
and actuators by the response regulator CheY-P (Figure 1A) affects 
the chemotactic performance of individual cell phenotypes. What 
CheY-P concentration maximizes drift velocity of cells navigating 
exponential gradients of methyl-aspartate? We examine this 
question using stochastic simulations of individual cells and an 
analytical model. 

Simulations were conducted using a standard model of the 
chemotaxis pathway in individual cells [2,15] as described in 
Methods. Receptor-kinase complex activity is modeled as an all- 
or-none response using quasi-equilibrium dynamics for fast ligand 
binding, chemoreceptor conformational changes, and phosphor- 
ylation cascade [15]. The slower (de)methylation kinetics follow 
simple negative integral feedback dynamics with adaptation rate 
T~ . The flagellar motor is modeled as an inhomogeneous Poisson 



process that switches cell behavior between runs and tumbles with 
rates defined as a function of CheY-P concentration that varies in 
time, T(t). The parameters of the motor model are calibrated to 
recent experimental measurements [19,20,23,24]. While motor 
adaptation [23] is not included at first, its effects are analyzed later 
in the paper. During runs, a cell swims with constant speed 
v = 20 u.m 1 in a direction subjected to rotational diffusion 
(rotation diffusion constant, D, = 0.062 rad 2 s 1 [1]). For simplic- 
ity the effects of multiple flagellar motors [2,35] or directional 
persistence [1] are not included but discussed in the Discussion 
section. Hence, in this model, motor clockwise bias (CW) and 
cell tumble bias (TB) are the same. We consider cells containing 
only Tar receptors and use methyl-aspartate as the ligand. 
Our results readily extend to more complex receptor cluster 
configurations. 

Three-dimensional trajectories of individual cells were simulat- 
ed as described in [2] for various cell phenotypes, which are 
characterized by the receptor adaptation times % and adapted 
CheY-P concentrations T 0 , in gradients of chemoattractant of 
different steepness, g. Following previous studies [25-27], we used 
exponential gradients (L(x) — Lg^") so that cells experience an 
approximately constant "force" from the attractant field, as the 
chemotaxis system is a fold-change detector (Eq. (2) below). This 
makes it possible to define a steady-state drift velocity, making the 
problem analytically tractable. The performance of each cell 
phenotype, which is defined by a unique adapted CheY-P 
concentration (T 0 ) and receptor adaptation time (t), in each 
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gradient steepness (g) is defined as the drift velocity Vd(Y 0 , t, g) 
along the gradient direction calculated by averaging the velocity of 
10,000 phenotypically identical cells over 4 minutes (Methods). 
The first simulations were done in a relatively shallow 
gradient g 1 = 5,000 Mm, with adaptation times of x = 5, 10, and 
30 s, and adapted CheY-P concentrations spanning the range 
T 0 = 1-4 uM. 

Plotting drift velocity as a function of the adapted CheY-P 
concentration reveals that maximal drift velocity is not achieved 
for CheY-P concentrations in the linear range of the CW bias 
response curve, where fluctuations in CheY-P result in large 
changes in clockwise bias (Figure 2A). Instead, it occurs when 
adapted CheY-P is at the lower end of this curve (around 2.4 M.M 
in Figure 2A). 

Analytical model of the drift velocity as a function of 
CheY-P concentration 

To understand the underlying reasons of this result, we derived 
an analytical relationship between CheY-P concentration and drift 
velocity along a one-dimensional gradient. For simplicity we used 
a one-dimensional analytical representation of bacterial chemo- 
taxis in two or three dimensions [27,33,34,36,37]. In this 
framework, cells either go up or down the gradient or tumble 
and the effect of rotational diffusion can be represented as a jump 
process between runs up and runs down the gradient with 
transition rate (d-l)D n where d represents the number of spatial 
dimension [36]. 

At quasi-steady state (for time scales longer than single run 
durations) and with no directional persistence (equal probability to 
run up or down the gradient), the drift velocity is proportional to 
the cell swimming speed v times the difference between the 
expected run durations up (t\ + iy R and down — 1} R the 
gradient divided by the total time including the time spent 
tumbling <?>r [36]: 



v <*l + l>*-<f|-l>* 
<f| + l + 4- 2<f> r 



(1) 



The only difference between d=2 or d=3 dimensions is a 
rescaling of the drift velocity and rotational diffusion (factor d in 
the equations above; see Text SI). 

The expected run duration up or down the gradient is 
controlled by the cellular concentration of CheY-P, T. This 
quantity is in turn a function of the free energy difference between 
the inactive and active receptor complexes, F, such that F=ln(a/ 
T-l), where a is the gain of the phosphorylation cascade. The 
receptor activity follows simple spring-like dynamics around the 
adapted free energy difference F 0 with adaptation time x 
(Methods, our Results still hold when considering asymmetric 
methylation/demethylation rates, see Text SI and Figure SI): 



dF 
dt 



(F-F Q )+sf 



(2) 



Here s = ± 1 or 0 for cells running up, down the gradient, or 
tumbling, respectively. As the cell moves along a trajectory x(t) it 
encounters different concentrations of the ligand L{x(t)). 
f = vN8 x \n[(l+L/K i )/(l+L/K a )] represents the magnitude of 
the change in free energy difference and depends on the local 
steepness of the gradient at the cell position. Here, v is the speed of 
the cell when it is swimming, JV is the gain of the cooperative 
receptor system and Ki and K a are the dissociation constants 
between ligand and receptors in the inactive and active 



conformation. In general, both s and f change as a function of 
time. For ligand concentrations Ki«L«K a we have 
/ « v N 8 X In L. If in addition, the gradient of ligand is exponential, 
L = L(j#" ' , we see that f~vNg becomes constant where g represents 
the inverse length scale of the gradient. Therefore, in an 
exponential gradient the free energy difference of the receptors, 
F, tends to increase at the constant positive rate +vNg when the cell 
swims up the gradient and to decrease at the constant negative rate 
— vNg when the cell swims down the gradient. This in turn causes 
the CheY-P concentration to decrease (increase) when cells swim 
up (down) the gradient. The exact CheY-P concentration 
trajectories can be calculated by integrating Eq. (2) as a function 
of time for different initial condition F t (Figure 2B), while a cell is 
swimming up (s = l green curve) or down (s=—l red curve) a 
gradient. 

The expected durations of a run, l ]i l (Y), or a tumble, lj '(JO, 
are plotted as a function of CheY-P concentration in Figure 2B 
(dashed lines; see definition in Methods). When a cell runs up or 
down the gradient, the rates of switching from one state to another 
change as a function of time and the direction of motion because 
they depend on the CheY-P trajectory. A run up or down the 
gradient can also be terminated by random reorientation from 
rotational diffusion with rate (d-\)D T [36]. Altogether, the rate at 
which a run is terminated by either rotational diffusion or a tumble 
is thus 1 = (d - 1 )D r + a r ( Y(F)). 

In a shallow gradient, F deviates little from the adapted value F 0 
and the adapted value of x R provides a good approximation of the 
expected run duration along a direction (black line in Figure 2B): 
x RO = x R (F 0 ) = ((d-l)D + 2. R o)~ 1 , where 1 R0 = 1 R (F 0 ). When 
swimming up or down the gradient, CheY-P fluctuates 
(Figure 2B, green lines for up, red lines for down) and the run 
lengths are modulated approximately following x R0 (black line and 
red and green circles in Figure 2B). According to equation (1), drift 
velocity is largest where the contrast between run duration up and 
down the gradient is the largest. Figure 2B reveals that this is the 
case where the slope of the expected run length as a function of 
CheY-P concentration is largest, which corresponds to the foot of 
the motor CW bias curve (Figure 2A) in agreement with the 
simulations. In contrast, for higher valued of CheY-P that are 
within the dynamic range of the CW bias response function, (e.g. 
T 0 = 3 uM in Figure 2B) run durations up and down the gradient 
have a smaller contrast and longer tumble duration (dashed line 
Figure 2B), resulting in slower drift velocity. 

In the limit of shallow gradients, Equation [1] can be linearized 
around the adapted values F 0 and Y 0 to obtain the drift velocity 
(Methods and Text SI): 



- t/%R0 



T«0 



fdtr 



x R0 ' (l-TB 0 )v 2 Ng 
1+t«o/t d 



(3) 



Here, TBq = a R o/(Xto + Aro) represents the tumble bias of the cell 
as a function of the adapted CheY-P concentration T 0 and the 
subscript 0 indicates that the rates X R and X T and x R = dx R /dF are 
all evaluated at the adapted state. The integral in Eq. (3) is the 
time-averaged input over the run durations, which in this 
approximation are exponentially distributed with characteristic 
time scale x R o- For Ki«L«K a and exponential gradients, the 
rate of change of the free energy difference sfxsvNg is constant 
during a run. Equation (3) indicates that the drift velocity is 
proportional to the gradient steepness g and the gain of the 
receptor cluster N. From Equation (3) we also obtain the 
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chemotaxis coefficient of an individual cell phenotype: 
X(Y 0 ,T)=V D (Y 0 ,T)/g. 

Plotting the drift velocity as a function of T 0 on top of the 
simulation results in Figure 2A shows that Equation (3) provides a 
good prediction of the drift velocity in shallow gradients and 
confirms that maximum velocity is reached for CheY-P values at 
the foot of the CW bias response curve. 

In this linear regime, the optimal CheY-P concentration is only 
weakly dependent on the cell adaptation time and does not depend 
on the gradient steepness. The factor Zrq /(I + T/m/t) in Equation 
(3) encapsulates the relationship between drift velocity and the 
CheY-P concentration (or free energy difference Fj. For small 
adaptation times, i«1ro, it increases linearly with adaptation 
time and is maximum where the slope dlriTRo/dFo is largest. 
For larger adaptation times, t»Tr 0 , this factor becomes 
~%rq = Yq( Yq/ol — \)dzjio/dY. Because the response function 
of the motor arq = ar(Yo) (defined in Methods) is very steep 
(dashed line in Figure 2B), the slope dtRo/dYo (slope of black line 
in Figure 2B) changes much faster as a function of Y 0 than Y 0 
(To/a-1). 

Because rotational diffusion imposes an upper bound on the run 
length along a given direction [27,32], it determines, along with 
the motor parameters, the optimal range for CheY-P fluctuations. 
As rotational diffusion becomes smaller, cells are able to maintain 
their original direction for a longer time. The upper bound on the 
run length therefore becomes longer (dashed lines in Figure 2C) 
and the optimal CheY-P concentration becomes smaller (full lines 
in Figure 2C). 

Changes in the switching frequency of the flagellar motor, CO, 
(see Methods) also affects the optimal CheY-P concentration. This 
becomes apparent when considering that the rate of switching 
from run to tumbles scales linearly with the switching rate of the 
flagellar motor, CO (see Methods). Therefore, the expected run 
duration of a cell scales like the inverse of Q). The result of this 
scaling is that for increasing values of CO the inflection point of the 
expected run length as a function of CheY-P shifts to lower values 
of CheY-P (dashed lines in Figure 2D). Thus, increasing the 
switching rate of the flagellar motors tends to decrease the optimal 
CheY-P concentration (full lines in Figure 2D). It also increases the 
maximum drift velocity that can be reached. 

The analytical model of drift velocity in Eq. (3) is different from 
previous results [27] in two ways. First, it takes into account both 
the adaptation time and the tumbling state of the cell. Taking the 
limit T— *oo,TBq xO in our model we recover the previous results. 
Second, in the previous study, the switching rate of the motor l R 
was a steep function of kinase activity centered at the adapted 
kinase activity level, or equivalently, at the adapted CheY-P 
concentration Y 0 . This means that changing Y 0 would also change 
the set point of the motor. However, the adapted CheY-P 
concentration and the set point of the motor are independent 
parameters. For this reason here we focused on the relationship 
between the set point of the motor and CheY-P activity. According 
to Eq. (3) the flagellar motor has its own sensitivity set point 
independent of the adapted CheY-P concentration of the sensory 
system Y 0 (see definition of A R in Methods; motor adaptation 
[23,38] is considered below). 

A unique operational CheY-P concentration maximizes 
drift velocity for multiple gradients and adaptation times 

Experiments have shown that when immobilized cells are 
exposed to an exponential ramp of methyl-aspartate, CheY-P 
activity reaches a new steady-state, Y m , which is lower than its 
adapted activity, Y 0 , because of the relatively slow adaptation rate 
of the system [9,26] (Figure IB second line). When cells are 
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Figure 3. Feedback of the behavior of cells swimming in 
exponential gradients onto the operational CheY-P concentra- 
tion. A. Temporal profiles of the average methyl-aspartate concentra- 
tion encountered by cells swimming in a steep exponential gradient 
(g _1 = 1,000 |rm). Different phenotypes are considered (solid black: 
y 0 = 2.4uJvl, t = 10s, solid gray: V 0 = 2.4frM, t = 30 s, dotted black: 
Y 0 = 3 u.M, x=10 s) (the y-axis is on a log scale). B. Corresponding 
average CheY-P concentration as a function of time in these same cells 
C. Magnitude of the drop in average CheY-P activity (difference 
between adapted and operational CheY-P concentrations Y m -Y 0 ) as a 
function of the drift velocity. Two different adaptation times are 
considered (black: x=10 s, grey: t = 30 s). The gradient is the same 
gradient as in panel A. Dots are averages over 10,000 stochastic 
simulations for populations with different adapted CheY-P concentra- 
tions (V 0 >2.4 JJ.IV1 in both cases). Lines are from Eq. (4). D. Drift velocity 
V D as a function of adapted CheY-P concentration, V 0 (filled circles), and 
operational CheY-P concentration, Y m (open circles) in stochastic 
simulations (average over 10,000 replicates for each circle, t = 10s). 
Y m is instantaneous CheY-P concentration averaged over the popula- 
tion while drifting between f=60 and 300 s). Two exponential 
gradients of methyl-aspartate are considered (g~ 1 = 1,000 |im (black), 
5,000 |im (grey)). Black arrow: cell population in blue in Figure 4C. 
doi:10.1371/journal.pcbi.1003694.g003 

swimming in an exponential gradient (Figure IB third line), we 
expected a similar effect to take place because the average drift of 
an individual cell up the gradient will cause this cell to experience, 
on average, an exponential increase in ligand concentration as it 
makes its way up the gradient. While this effect should be minimal 
in a shallow gradient, it could become important in steep or 
rapidly changing gradients [27,39], especially for cells with longer 
adaptation times. 

To investigate this issue we simulated cells swimming in a 
steeper exponential gradient (g _1 = 1,000 urn). After less than one 
minute of simulation, cell populations (10,000 replicate trajectories 
for each phenotype) reached a constant steady state drift velocity. 
We calculated the average ligand concentration that the cells 
encountered over time (Figure 3A). This reveals that the swimming 
cells experience an average exponential increase in ligand 
concentration over time. This average input is similar to the 
signal experienced by immobilized cells exposed to temporal 
exponential ramps (Figure IB, second line) [9,26]. However, for 
the swimming cells the ramp rate is dynamically determined by the 
average drift velocity in the direction of the gradient (Figure 3A). 
Thus, for swimming cells the ramp rate depends on the feedback 
of the performance onto the input signal (Figure IB, third line). 
Consistent with experimental results obtained with immobilized 
cells exposed to exponential ramps [26], the average CheY-P 
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concentration in the swimming cells reaches a stable dynamical 
equilibrium, the operational value Y m , after an initial drop from 
the adapted CheY-P concentration (Y 0 ) (Figure 3B). 

The fact that the operational CheY-P concentration is not the 
same as the adapted CheY-P concentration implies that an 
optimal choice of adapted CheY-P must take into account this 
behavioral feedback. While a phenotype may for example have an 
adapted CheY-P concentration equal to the optimal concentration 
(~2.4 |J.M), during chemotaxis this level drops to an operational 
level lower than the optimum, hindering its performance 
(Figure 3B, solid black line). This effect is intensified when the 
adaptation time of the receptor cluster increases (Figure 3B, grey 
line). On the other hand, a phenotype with an adapted CheY-P 
concentration higher than the optimal concentration can ap- 
proach the optimal operational CheY-P concentration as it 
reaches its steady-state drift velocity (Figure 3B, dotted line). The 
difference between Y 0 and Y m grows larger as drift velocity or the 
receptor adaptation time increase (Figure 3C). 

To determine whether the adapted or the operational CheY- 
P concentration is the primary variable that controls the 
average drift velocity in exponential gradients, we simulated 
cell populations with different adapted CheY-P concentrations 
and calculated their respective operational CheY-P concentra- 
tions. In a steep gradient (g -1 = 1000 |im), the optimal adapted 
CheY-P increased to ~2.7 |J,M compared to ~2.4 uM in a 
shallow gradient (Figure 3D). However, the optimal operational 
CheY-P concentrations for steep and shallow gradients are 
identical (Figure 3D). This suggests that a unique operational 
CheY-P concentration maximizes drift velocity in multiple 
gradients. 

The situation is different when the feedback is strong. In this 
case the signaling pathway fluctuates around the operational 
values F m and Y m , rather than the adapted values F 0 and Y 0 , 
Therefore, we need to update the analytical model to describe the 
drift velocity, V D , as a function of F m . If we linearize the drift 
velocity equation around F m rather than F 0 we obtain an equation 
identical to Eq. (3) but with the subscript 0 replaced by m and iR m , 
l'ltm> ^Rm, ^Tm, an d TB m now functions of F m . Knowing how V D 
depends on F m is not enough to calculate the drift velocity. We also 
need an equation that describes how F„ depends on Vo- To model 
the effect of a constant drift velocity along the chemical gradient 
on the activity of the receptor cluster, we can expand Eq. (2) 
around F m and solve for quasi-steady state: 



-- Fq + V D -cN-r In 
ax 



l+L(x(i))/Ki\ 

Y^Mkr Fa+VDXNs (4) 



This expression quantifies the deviation between the operational 
free energy difference F m and the adapted free energy difference F 0 
as a function of the drift velocity and is consistent with the results 
of our simulations (Figure 3C) and [27]. Equation (4) also makes 
clear that the strength of the feedback depends on adaptation time, 
the receptor cluster gain jV, and the steepness of the gradient. 
Behavioral feedback strongly affects performance because it moves 
Y m away from the optimal operating point relative to the motor. 
This, in turn, affects the capability of the motors to best use the 
information carried by CheY-P. 

By explicitly taking into account the effect of the behavioral 
feedback onto the coupling between the operating regime of 
CheY-P and the motor, Eqs. (3) (with g_^„) and (4) extend previous 
studies [17,25,27,32-34,36] and reveal new possible dynamical 
regimes for the biased random strategy as shown below. 
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Figure 4. Behavioral feedback can create a chemotactic "trap". 

A. Analytical drift velocity, V D as a function of Y m (Eq. (3) with o-^m, solid 
line) and feedback of V D on Y m (Eq. (4), dashed line) (t = 5s, 
g _1 = 1,000 urn, Y 0 = 2.7 \iM). Steady state drift velocity {Y m = 2A8 uM, 
black circle). B. Same as panel A, but for cells with longer adaptation 
time and higher adapted CheY-P (t = 30 s, V 0 = 3.5 uM). V d has three 
possible steady states: two stable (Y m2 = 2A aM and V m , = 3.49 llM 
(black dots)), and one unstable (V m = 2.97 uM, white dot). C. Individual 
drift velocities (in the direction of the gradient) and root mean square 
displacements (perpendicular to the gradient) of two different 
populations of 10,000 simulated cells (blue: r = 10 s, V 0 = 2.6 uM, red: 
t = 30 s, Y 0 = 3 |jJVI). D. Average V D as a function of Y 0 (filled circles) and 
Y m (open circles) for cells with a long adaptation time (t = 30 s). Black 
arrow: cell population plotted in panel C (red). 
doi:10.1371/journal.pcbi.1003694.g004 



Strong behavioral feedback can push the system through 
a bifurcation creating two possible chemotactic states for 
some cell phenotype: A fast drift state and a trapped 
state 

For a given phenotype (To, t) and gradient length-scale, the 
steady state drift velocity is determined by the intersection of two 
curves (Figure 4A). The first curve (solid line in Figure 4A) 
describes how the drift velocity depends on the operational CheY- 
P concentration, Y m . It is defined by Equation (3) (with fJ ^ m ) and its 
profile can be interpreted as follows. For very low values of CheY- 
P the cell never tumbles. Thus, the cell diffuses equally in all 
directions and the net drift along the gradient is zero. For high 
values of CheY-P, the cell tumbles all the time so drift is zero as 
well. In between these two extremes, drift velocity is maximized for 
a specific value of the operational CheY-P concentration. 
However, Y m is not an independent variable. As we showed above 
(Eq. (4)), because of the feedback the behavior onto the input, 
the operational CheY-P concentration is itself a function 
of the drift velocity (which can also be written as: 
Y m = a/(\+e xN s VD (a./Y 0 -X))). This equation defines the 
dashed line in Figure 4A, which intersects the horizontal axis at 
Yg. Because each line in Figure 4A defines a relationship between 
Vo and Y m , the intersection between the two lines fully determines 
the drift velocity and the operational CheY-P concentration (black 
circle in Figure 4A) for a given phenotype and gradient. 
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When the feedback is weak (zNg small, i.e. short adaptation 
time, small gain, or shallow gradient), the operational CheY-P 
concentration only exhibits a weak dependency on drift velocity 
and there is only one steady-state solution (intersection). There- 
fore, an appropriate adapted CheY-P concentration could be 
selected to ensure that operational CheY-P concentration is 
approximately optimal at all times (Figure 4A). 

When the feedback is stronger, drift velocity always acts as a 
negative feedback onto the operational CheY-P concentration. In 
contrast, the effect of the operational CheY-P concentration onto 
drift velocity depends on whether the operational CheY-P 
concentration is below or above the CheY-P concentration that 
maximizes drift velocity (Figure 4A). Below this concentration, the 
system obeys negative feedback dynamics, whereas above it, the 
system obeys positive feedback dynamics. This positive feedback 
loops combined with the non-linear decrease of the drift velocity as 
a function of the operational CheY-P concentration, which arise 
from the extreme sensitivity of the flagellar motor, can lead to 
bistability [40] . Indeed, for a stronger feedback (steeper gradient or 
longer adaptation time) the slope of the feedback curve (dashed 
line in Figure 4AB), which is proportional to 1 / rjVjj, decreases. 
Thus, for phenotypes with high enough adapted CheY-P 
concentration (T 0 is the intersection of the dashed line with the 
horizontal axis), the two curves can intersect more than once 
(Figure 4B). In this case, a single phenotype can now experience 
three different chemotactic states. Two of these states, Y mI and Y m2 
(filled circles in Figure 4B), are stable and are separated by one 
unstable state (open circle in Figure 4B). For one of the stable 
solution, the drift velocity is nearly zero and Y ml is high and very 
close to the adapted CheY-P concentration. For the other stable 
solution, the drift velocity is large and Y m 2 is much smaller than the 
adapted CheY-P concentration. 

This analysis suggests that an individual phenotype might 
experience two different chemotactic states with dramatically 
different performance: a fast drifting state and a "trapped" state. 
To find evidence of these two behaviors, we simulated two cell 
phenotypes (10,000 replicates for each phenotype) in a steep 
exponential gradient (g 1 = 1,000 urn). One phenotype was 
predicted to operate closer to the bifurcation than the other (red 
and blue dots in Figure 4C, respectively). Although both 
phenotypes reached the same average operational CheY-P 
(Y m — 2.3 uM), cells with a phenotype closer to the predicted 



bifurcation point (Y 0 = 3 |0,M, T = 30 s) exhibited a distribution of 
behavior (both drift velocity and diffusion) significantly skewed 
toward the "trapped" state (Figure 4C). Closer examination of the 
trajectories and CheY-P dynamics of individual cells in this 
simulation reveals that individual cells transition stochastically 
back and forth between the "trapped" and fast drifting state 
(Figure S2). For cases with higher feedback strength cells spend 
more and more time within the "trapped" state. 

When the feedback is strong and the system becomes multi- 
stable, the average includes cells in both the "trapped" and high 
drift states. This phenomenon explains the decreased average drift 
velocities observed when adaptation time is increased (above 
10 seconds) in a relatively steep gradient (g~' = 1,000 |J.m). It also 
explains the resulting shift of the best operational CheY-P 
concentration to lower concentrations (from ~2.4 to 2.1 \iM in 
Figure 4D), since for phenotypes with lower values of the adapted 
CheY-P only one stable state exists. Similar results are obtained 
when asymmetric methylation/ demethylation rates of the recep- 
tors are taken into account (Figure S3). 

Motor adaptation partially alleviates the chemotactic 
"trap" 

Recent experiments have shown that the number of FliM 
monomers in the C-ring of the flagellar motor slowly (~minutes) 
adapts as a function of the CW bias, affecting both the steepness 
and the half-maximum CheY-P concentration of the CW bias 
motor response curve [24]. To examine the effect of motor 
adaptation on the relationship between CheY-P concentrations 
and drift velocity, we added motor adaptation to our stochastic 
model of an individual chemotactic cell by taking into account 
recent experimental data [20,23,24] (Methods). The resulting CW 
bias response curve of the adapted motor agrees well with both 
recent [20] and earlier [19] experimental measurements. In fact, it 
matches earlier experiments [19] better than a simple Hill 
function, suggesting that in these experiments the individual 
motors measured had adapted to the particular concentration of 
CheY-P expressed in the corresponding individual cells (Figure 5A; 
Methods). 

Simulations of cells with motor adaptation in a shallow gradient 
{g 1 = 5,000 u.m) show that motor adaptation changes the shape of 
the drift velocity curve as a function of operational CheY-P, 
especially at high CheY-P concentrations (compare Figures 5B and 




Figure 5. Effect of motor adaptation on drift velocity V D in exponential gradients. A. Motor CW bias response curve as function of CheY-P 
concentration when the motor is allowed to adapt (solid line) fitted to data from [19] (circles; derivation in Materials and Methods). B. Average drift 
velocity as a function of operational CheY-P concentration Y m , in a shallow gradient. Same adaptation times and gradient steepness as Figure 2A. 
Lines: analytical solutions; circles: stochastic simulations (averages between f= 10 and 15 min are used to calculate V D (V m )). C. Same as Figure 4B but 
with motor adaptation. The drift velocity has only one stable steady sate (V m = 1.6 uM, black dot). Motor adaptation eliminated the other states 
present in Fig. 4B. 

doi:10.1371/journal.pcbi.1003694.g005 
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Adaptation time (s) Adaptation time (s) 

Figure 6. Performance trade-off in bacterial chemotaxis. A. Optimal adapted CheY-P concentrations Yq (solution of Eqs (4) with o_^. m and (5)) 
as a function of the chemoreceptor adaptation time in different exponential gradients (g _1 = 1,000 (red), 2,000 (green), and 5,000 (blue) urn). Dots 
indicate when the maximal theoretical drift velocities cross the bifurcation point (dotted lines represent the inaccessible optimal state). The optimal 
operational CheY-P concentration Y m is identical for all gradient length scales (black dashed line). B. Contour plot of drift velocities as a function of 
adaptation time and the adapted cell tumble bias in different exponential gradients (same colors as A). 75%, 90%, and 95% contours of the maximal 
theoretical drift velocities for each gradient (colors intensities from light to dark). Black dot: the best cell phenotype that achieves equal relative drift 
velocities in all three gradients (60% of the maximal V D with t = 7.5 s and TB 0 = 0.044). 
doi:10.1371/journal.pcbi.1003694.g006 



2A). These results are predicted by the analytical model (Eq. (3) 
with 0 _^ m ) once modified to include motor adaptation (Methods; 
lines in Figure 5B). Setting the adapted activity of the motor for a 
given CheY-P concentration to lower or higher CW biases gives 
qualitatively equivalent results (Figure S4 and S5). 

How does the motor adaptation affect the bifurcation? In a 
steep gradient the behavioral feedback (Eq. (4)) must be taken into 
account (Figure 5C dashed line). Comparing Figure 5C and 4B we 
see that motor adaptation enable cells with high adapted CheY-P 
concentration to avoid the chemotactic trap improving perfor- 
mance (see Figure S6). This should provide a selective advantage 
because it helps buffer the functional consequences of inevitable 
cell-to-cell variability in the adapted CheY-P concentration, by 
increasing the range of CheY expression levels that allows effective 
chemotaxis. 

Motor adaptation also affects the optimal operational CheY-P 
concentration (compare Figures 5B and 2A), shifting it to lower 
concentrations. When cells are drifting up a gradient, CheY-P 
drops to the operational CheY-P, causing the CW bias to drop. 
With motor adaptation, the lower operating CW bias causes the 
motor to shift its sensitivity to a lower CheY-P concentration. We 
see again that maximal performance is reached for T„ at the 
bottom of the CW bias response curve of the now adapted motor 
(Figure 5A). However, the motor can only compensate partially for 
the shift in operational CheY-P concentration. 

An individual phenotype faces a performance trade-off in 
different gradients 

As long as the system does not undergo bifurcation, maximum 
drift velocity is achieved by having a long adaptation time while 
maintaining the operational concentration of CheY-P in the 
optimal range. Therefore, the optimal adapted CheY-P concen- 
tration depends on the gradient length-scale and the adaptation 
time (Figure 6). 

In shallow gradients, the strength of the feedback is small, as is 
the difference between operational and adapted CheY-P. Thus, it 
is possible to select an adapted CheY-P concentration that will 
perform relatively well for multiple adaptation times (Figure 6A 
blue line). In steeper gradients, the feedback is stronger (Eq. (4)) 
and the difference between T m and Y 0 grows larger with adaptation 
time. Maintaining the optimal operational CheY-P concentration 



requires a higher adapted CheY-P concentration (Figure 5A green 
and red). The bifurcation of the system imposes an upper bound 
on the range of Tg beyond which a portion of the cells spend a 
significant amount of time trapped into a non-optimal state even 
with motor adaptation (Figure 6A dashed lines). Therefore, the 
optimal adapted CheY-P concentration is a function of both 
receptor adaptation time and gradient length-scale, making it 
difficult for a single phenotype to maximize drift velocity in 
multiple environments (Figure 6A). 

To characterize the resulting performance trade-off and map it 
to phenotypic space, we calculated the contours of drift velocity 
relative to its maximum in each environment, as a function of 
adaptation time T and the adapted cell tumble bias (Figure 6B). In 
shallow gradients, cells benefit from a relatively long adaptation 
time and a low adapted CW bias. In steep gradients, cells benefit 
from a short adaptation time and a higher adapted CW bias. The 
best generalist phenotype can achieve at most 60% relative 
performance in all three gradients considered here. Motor 
adaptation, which was taken into account in generating Figure 6, 
alleviates only partially the tradeoff faced by single cells. 

Discussion 

The adaptive response and feedback control of the receptor 
cluster play a critical role in the robustness of the chemotaxis 
system [8,10,15,17,33]. However, chemotactic performance also 
relies on the optimal operation of the flagellar motors, which 
directly control cell behavior. By focusing on how the CheY-P 
concentration affects the coupling between sensors and actuators, 
we revealed the existence of an operational regime for CheY-P 
concentration, which is distinct from the adapted CheY-P 
concentration, that maximizes drift velocity in a wide range of 
gradient length-scales and receptor adaptation times. Fluctuations 
around the best operational CheY-P concentrations maximize the 
contrast between run duration up and down the gradients. This 
occurs outside the most sensitive region of the CW bias response 
curve of the motor. Thus, chemotactic performance relies on 
maintaining the operational CheY-P concentration within bounds 
[21,22,41] around this optimal value. 

The best operational CheY-P concentration is also determined 
by the cell rotational diffusion constant D n which imposes an 
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upper bound on run durations in any particular direction 
(Figure 2C) [27,32]. In a more viscous environment or for longer 
cells, the lower rotational diffusion will result in a lower optimal 
operational CheY-P concentration. For an ellipsoid, rotational 
diffusion is inversely proportional the length of the major axis. 
Therefore, as cells grow the optimal range will shift to lower 
CheY-P concentrations. If the cell maintains a constant amount of 
CheY as it grew, the effective concentration of CheY-P would 
decrease, resulting in robust chemotactic performance during cell 
growth. 

The switching frequency of the flagellar motors also affects the 
best operational CheY-P concentration. Higher switching fre- 
quencies tend to increase drift velocity while shifting the maximum 
to smaller CheY-P concentrations (Figure 2D). Therefore, the best 
operational CheY-P concentration is further away from the motor 
threshold. However, the range of CheY-P concentrations where 
the drift velocity is high becomes narrower (because the expected 
run length becomes a steeper function of CheY-P). This tends to 
increase the performance trade-off between different gradient 
length-scales. Thus, while selecting a higher switching frequency 
for the flagellar motors may improve performance of some 
phenotypes it may be detrimental for the population overall. 
Another important consideration is that the switching frequency is 
bounded by the speed at which the motor and associated flagella 
can switch confirmation [2,42]. 

Directional persistence (amount by which the swimming 
direction of a new run is biased towards the swimming direction 
of previous run) has been shown to affect chemotactic perfor- 
mance in climbing shallow gradients of attractants [1,43-45]. 
However, previous modeling and simulations efforts have been 
done using cells with non-optimal CheY-P concentrations (usually 
at 3 U.M). In this regime, cells have a high tumbling rate, short run 
lengths, and low drift velocity. Directional persistence effectively 
reduces the reorientation rate of cells [45], which is equivalent to 
reducing the tumbling rate slightly. Therefore, directional 
persistence will shift the optimal CheY-P concentration to higher 
concentrations and improve the drift velocity of frequendy 
tumbling cells [44]. On the other hand, when cells operate at or 
close to the optimal CheY-P concentration, the tumbling rate is 
low. Therefore, the run length in a given direction is terminated by 
rotational diffusion and not by tumbles. For optimal phenotypes, 
the relative effect of directional persistence on chemotactic drift is 
thus less important. 

Previous studies have examined how the adaptation time affects 
chemotactic performance [12,25,27,32,46]. However, these stud- 
ies only considered single values for the adapted CheY-P 
concentration (typically set to a CW bias of 0.5) and concluded 
that adaptation time should decrease as gradients get steeper to 
keep the operational CheY-P concentration within the dynamic 
range of the motor CW bias response. We found that, as long as 
the cell can maintain the optimal operational CheY-P concentra- 
tion, longer adaptation time is better because it enhances input 
signal over the course of a run. However, long adaptation 
reinforces the feedback from the cell drift velocity on the system 
and can lead to undesirable bistability. Therefore, the bifurcation 
boundary imposes an upper limit on adaptation time as a function 
of the gradient length-scale. Interestingly, the distribution of 
tumble bias typically observed during exponential growth in single 
E. coli cells ranges from 0. 1 to about 0.4 and not many cells are 
found that have higher tumble bias [31]. Selection for cells with 
tumble bias below 0.4 is consistent with our finding that the 
performance of cells with higher tumble bias will suffer from the 
existence of the "trapped" chemotaxis state. 



Our results also provide a strong justification for the role of the 
recently-discovered flagellar motor adaptation. Indeed, we found 
that motor adaptation [23,38] plays a significant role in mitigating 
the behavioral feedback for cells with high tumble bias. When such 
feedback was included, cells with high tumble bias could escape 
the "trap" and gain access to a high drifting state in steep gradient. 
Our model also resolved an apparent contradiction between two 
sets of experimental measurements of the CW bias response of the 
flagellar motor as a function of CheY-P concentration. While one 
measurements reported a Hill coefficient of n = 10 [19], newer 
experiments reported a Hill coefficient of n = 20 [20]. In this paper 
we used the new value n = 20 and showed that the previous 
measurements are fitted with the same parameter value if one 
makes the reasonable assumption that the motors had had time to 
adapt before each individual cell measurement (Figure 5A). 

Because the difference between the operational and adapted 
CheY-P concentrations depends on the strength of the behavioral 
feedback, which itself is proportional to gradient steepness, 
different adapted CheY-P concentrations and adaptation times 
are required to perform optimally in different gradients. Thus, in 
conditions where drift velocity is important, cells are faced with a 
performance trade-off. Even though motor adaptation was 
included, the best compromising phenotype over the gradient 
steepness considered in this study achieved at most 60% of the 
theoretical maximal drift velocity in all gradients. The observed 
cell-to-cell phenotypic diversity in adaptation time and adapted 
tumble bias [29,31] in an isogenic population may resolve the 
performance trade-off faced by single cells to improve the chance 
of survival of a unique genotype in complex or varying 
environments. In addition, the negative correlation between 
tumble bias and adaptation time observed by Park et al. in an 
isogenic population of E. coli [31], is consistent with our 
predictions about the most beneficial way to distribute phenotypes 
(Figure 6B). 

At its core, the biased random walk relies on the dynamical 
control of the probability of reorientation. Overall, our analysis 
reveals limits to the use of negative integral feedback to control 
such strategy. Because the biased random walk strategy is used by 
many organisms, these results will inform our understanding of the 
constraints faced by other organisms as well. 

Materials and Methods 

Model and simulations 

We used a standard model of bacterial chemotaxis [15] as 
described in [2]. For a cell following the trajectory x(t), the output 
of the sensory module, the CheY-P concentration, is 
Y(t) = a/ ( 1 +e F ) where the free energy difference 
between inactive and active receptor complexes, 
E = £o + £im + N\n({\+L/K i )/(l+L/K a j), is a function of the 
methylation level, m(t) and ligand concentration L(x(t)). With 
a = 6 uM, e 0 = 6, s, = -1,N = 6, andA; = 0.0182 mM,K a = 3 mM 
for methyl-aspartate and Tar receptors in the inactive and active 
conformation. When the cell is adapted to its environment, 
Fq = £o +£iWo = ln(a/ To — !)■ Adaptation mediated by methyla- 
tion and demethylation of the receptors follows 
dm/dt= —(m — m(L))/x, where m — m{L)={F — Fq)/e\. The 
methylation level m is positive and bounded by the total 
number of methylation sites m max = 48 available in a cooperative 
unit of receptors. The resulting adaptation dynamics fits 
recent experiments [26]. Cells switch between runs, R, and 
tumble, T, with ratesA^T = (oexp[+G(Y(F))]. The motor 
is modeled as a bistable system with switching frequency 
eo = 1.3 s _1 (unless otherwise stated) and free energy difference 
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G(Y)=e 2 /4-(ei/2)(l+K/Y)~ 1 where s 2 and 8 3 are non- 
dimensional constants that control the basal rate of switching of 
the motor when Y= 0 and the degree of cooperativity of the 
motor, respectively. K is the binding constant of CheY-P to FliM at 
the base of the motor. With 8^ = 83 = 80, and K= 3.06 uM, this 
coarse-grained motor model fits well recent experimental mea- 
surements of CW bias (Hill coefficient 20) and switching frequency 
[20,23,24]. Motor adaptation is considered below. 

Linear expansion 

Eq. (2) follows by taking the time derivative of F and using the 
relations from the previous section. Integration of Eq. (2) gives: 
F{t,s,Fj) =F 0 + (F i -F 0 )e-'' t + se-<^ < e u lj(u)du. The expected 
duration of a run along the direction s = ±1 is determined by the 
integral of the rate x R 1 of terminating a run along the direction s 
by tumbling or because of rotational diffusion: 



dt 



(5) 



Because the average cell drift velocity in the direction of the gradient 
is determined by the contrast between expected run durations up and 
down the gradient (Equation (1)), the quantity of interest to calculate 
from Equation (5) is (t\ + \,Fi) — (t\ — l,F,y. In a shallow gradient, 
the deviations from the adapted free energy difference F 0 are small. 
Considering only first order deviations AF = F — Fq around F 0 the 
change in free energy \AF/Fq\bs a response to changes in ligand 
concentration is small and the inverse of the rate of run termination can 
be approximated by t r (F)xt R o + x' R qAF where the mean run 
duration along a direction Trq = t r (Y(F( ) )) = ((d— l)D r + Aro) 1 
and the gainr^o = dx R jdF are evaluated at F 0 . Similar linear 
expansions are carried out for 1 R and It- Linear expansion of 
the free energy difference in Eq. (5) and integration by part 
gives 0\s,F lR y ^r R0 [l - X R0 £ e-'^MAFlitls^dt] +0{AF 2 ). 

CO CO . 

For tumble, (t\F iT }» f <?~ A ™' f A R0 e-'R/ T R0dt R dt + O(AF) = 
0 b 

tRO^RO /^T0 + O( AF) . Inserting in Eq. (1) and using the solution F[t, s, 
Fj) we obtain the drift velocity to first order in AF (Eq. (3)). 

Motor adaptation 

The number of FliM molecules in the motor, n, is modeled as a 
binding and unbinding process with CW bias dependent 
rates [38]: dn/dt = k on {\ - CW) /{I + An/(n 2 -n)) -k off CW / 
(1 +An/(n — «i)).The constants k„„ and k 0 jr define the rate of 
adaptation of the motor. n\ and « 2 are the minimum and 
maximum FliM ring size that a motor can accommodate. An is an 
effective half max parameter that guarantees that the effective 
rates of unbinding and binding to the motor go to zero when n 
approaches ri\ or n 2 . When n changes it affects the steepness of the 
motor CW bias response, CW = 1/(1 +e 2G ), which in our case is 
controlled by E3 (see above). We used a simple linear relationship 
£3 = £3,1 (n — «o) +£3,0 where 83^ is the slope and n Q and 83^ are 
the pre-stimuli level of the number of FliM and motor steepness, 
respectively. k m = 0.025 s -1 , «i = 34, n 2 = 44 from [24]. We 
choose £3 0 = 80, no = 36 to match the Hill coefficient of 20 
measured for individual motor response curves [20], and fit 
An = 4. 16, 8 3j i = 1.96 to reproduce [19] (Figure 5B inset). 
k„jf= 0.0063 s 1 controls the CW bias that the motor adapts to 
(0.2 in this case, typical for wild type population of E. coli selected 
for swimming on agar plates [31]). At steady state, dn/dt=0 
defines CW{n(e 3 )) (Eq. [S20] in Text SI). On the other hand, 
assuming quasi-equilibrium between the motor and operational 



CheY-P concentration Y m , we have CW( Y m ,e 3 ) =1/(1 +e 2G( Y "" 63) ). 
Solving the two equations gives £3 as a function of Y m from which we 
can calculate the drift velocity as Eq. (4) with motor adaptation 
(Figure 5). 

Supporting Information 

Figure SI Effect of asymmetric methylation/demethylation 
rates on drift velocity Fd in exponential gradient. Simulated drift 
velocity Vo (average velocity of 10,000 cells between <=60 and 
300 s) as a function of operational CheY-P concentration Y m in a 
shallow gradient {L„- 200 \>M and g' 1 = 5,000 |lm) for cells with 
methylation rates V R = 0A s _1 (red), 0.2 s -1 (green), and 0.4 s -1 
(blue). 
(PDF) 

Figure S2 A simulated cell can transition in and out of the non- 

chemotactic state to reach the high drift velocity state 
when swimming in a steep gradient of methyl-aspartate 
(g _1 = 1,000 urn) illustrating the bi-stable behavior of this cell 
phenotype (Y 0 = 3.0 u.M and 1 — 30 s). A. Single cell drift velocity 
as a function of its operational CheY-P concentration. When the 
cell escapes the "trapped" chemotactic state, characterized by a 
high CheY-P concentration, the behavioral feedback maintains an 
optimal CheY-P concentration and a high drift velocity. B. Cell 
position along the gradient as a function of its operational CheY-P 
concentration. The cell can escape the low drift velocity state and 
maintain a low CheY-P concentration when running up the 
gradient. On the other hand, the cell can return to the "trapped" 
state after a long run down the gradient. The CheY-P 
concentration and drift velocity were calculated over a 
moving average window of 10 seconds. The time progression 
along the trajectory is indicated by the color of the stroke from 
blue to red. 
(PDF) 

Figure S3 Effect of asymmetric methylation/demethylation 
rates on drift velocity Vo in steep exponential gradient. Vp from 
stochastic simulations (methylation rate Vr = 0. 1 s ') as a function 
of Y 0 (filled circles) and Y m (open circles) in exponential gradient of 
methyl-aspartate (g 1 = 1,000 \im). 
(PDF) 

Figure S4 Drift velocity as a function of operational CheY-P 
when the rate of binding between FliM and the motor is 
k 0 ff= 0.025 s . Circles are from simulations. Lines are from 
analytical solution. Everything is the same as in Figure 5A. The 
only difference is the value of k 0 jp 
(PDF) 

Figure S5 Drift velocity as a function of operational CheY-P 
when the rate of binding between FliM and the motor is 
k 0 ff= 0.001 3 s In this case, the motor does not adapt fast 
enough to reach quasi-steady state. The analytical solution (lines) 
makes the approximation that the system is at steady state. 
(PDF) 

Figure S6 Scatter plot of individual drift velocities (in the 
direction of the gradient) and root mean square displacements 
(perpendicular to the gradient) of 10,000 simulated cells with 
motor adaptation (black) and without motor adaptation (red), with 
adaptation time t = 30 s adapted CheY-P concentration 
Y 0 = 3.5 |iM, gradient length scale g -1 = 1,000 um. 
(PDF) 



Text SI 

(PDF) 



Detailed derivation of the analytical model. 



PLOS Computational Biology | www.ploscompbiol.org 



10 



June 2014 | Volume 10 | Issue 6 | e1 003694 



Limits of Bacterial Chemotaxis 



Acknowledgments 

We thank D. Clark, N. Frankel, J Long, N. Olsman, A. Waite, W. Pontius, 
and S. Zuckcr for helpful discussions. Simulations were performed at the 
Yale University High Performance Computing Center. 



Author Contributions 

Conceived and designed the experiments: YSD XF LHN TE. Performed 
the experiments: YSD XF LHN. Analyzed the data: YSD XF LHN TE. 
Wrote the paper: YSD XF TE. 



References 

1. Berg HC, Brown DA (1972) Chemotaxis in Escherichia eoli analysed by three- 
dimensional tracking. Nature 239: 500—504. 

2. Sneddon MW, Pontius W, Emonct T (2012) Stochastic coordination of multiple 
actuators reduces latency and improves chemotaetic response in bacteria. Proe 
Natl Acad Sci U S A 109: 805-810. 

3. Mm TL, Mcars PJ, Chubiz LM, Rao CV, Golding I, et al. (2009) High- 
resolution, long-term characterization of bacterial motility using optical 
tweezers. Nat Methods 6: 831-835. 

4. Sourjik V, Wingrccn NS (2012) Responding to chemical gradients: bacterial 
chemotaxis. Current opinion in cell biology 24: 262—268. 

5. Bray E), Levin MD, Morton-Firth CJ (1998) Receptor clustering as a cellular 
mechanism to control sensitivity. Nature 393: 85—88. 

6. Sourjik V, Berg HC (2004) Functional interactions between receptors in 
bacterial chemotaxis. Nature 428: 437-441. 

7. Li M, Khursigara CM, Subramaniam S, Hazelbauer GL (2011) Chemotaxis 
kinase Che A is activated by three neighbouring chcmorcccptor dimers as 
effectively as by receptor clusters. Mol Microbiol 79: 677—685. 

8. Alon U, Surette MG, Barkai N, Leibler S (1999) Robustness in bacterial 
chemotaxis. Nature 397: 168-171. 

9. Block SM, Scgall JE, Berg HC (1983) Adaptation kinetics in bacterial 
chemotaxis. J Bacteriol 154: 312-323. 

10. Yi TM, Huang Y, Simon MI, Doyle J (2000) Robust perfect adaptation in 
bacterial chemotaxis through integral feedback control. Proe Natl Acad Sei USA 
97: 4649-4653. 

11. Korobkova E, Emonet T, Vilar JM, Shimizu TS, Cluzel P (2004) From 
molecular noise to behavioural variability in a single bacterium. Nature 428: 
574-578. 

12. Pontius W, Sneddon MW, Emonct T (2013) Adaptation dynamics in densely 
clustered chemo receptors. PLoS Comput Biol 9: el003230. 

13. Mm TL, Mears PJ, Golding I, Chemla YR (2012) Chemotaetic adaptation 
kinetics of individual Escherichia coli cells. Proe Natl Acad Sci USA 109: 9869- 
9874. 

14. Lazova MD, Ahmed T, Bellomo D, Stoekcr R, Shimizu TS (2011) Response 
rescaling in bacterial chemotaxis. Proe Nad Acad Sci USA 108: 13870-13875. 

15. Tu Y (2013) Quantitative modeling of bacterial chemotaxis: signal amplification 
and accurate adaptation. Annu Rev Biophys 42: 337—359. 

16. Neumann S, Vladimirov N, Krembel AK, Wingrccn NS, Sourjik V (2014) 
Imprecision of adaptation in Escherichia coli chemotaxis. PLoS One 9: c84904. 

17. Celani A, Vergassola M (2010) Bacterial strategies for chemotaxis response. Proe 
Natl Acad Sci USA 107: 1391-1396. 

18. Berg HC (2003) The rotary motor of bacterial flagclla. Annual review of 
biochemistry 72: 19-54. 

19. Cluzel P, Surette M, Leibler S (2000) An ultrasensitive bacterial motor revealed 
by monitoring signaling proteins in single cells. Science (New York, NY) 287: 
1652-1655. 

20. YuanJ, Berg HC (2013) Ultrascnsitivity of an adaptive bacterial motor. Journal 
of molecular biology 425: 1760—1764. 

21. Kollmann M, Lovdok L, Bartholome K, Timmer J, Sourjik V (2005) Design 
principles of a bacterial signalling network. Nature 438: 504-507. 

22. Lovdok L, Bentele K, Vladimirov N, Muller A, Pop FS, et al. (2009) Role of 
translational coupling in robustness of bacterial chemotaxis pathway. PLoS Biol 
7: el000171. 

23. YuanJ, Branch RW, Hosu BG, Berg HC (2012) Adaptation at the output of the 
chemotaxis signalling pathway. Nature 484: 233—236. 

24. Lele PP, Branch RW, Nathan VS, Berg HC (2012) Mechanism for adaptive 
remodeling of the bacterial flagellar switch. Proe Natl Acad Sei USA 109: 
20018-20022. 



25. Vladimirov N, Lovdok L, Lebiedz D, Sourjik V (2008) Dependence of bacterial 
chemotaxis on gradient shape and adaptation rate. PLoS Comput Biol 4: 
c 1000242. 

26. Shimizu TS, Tu Y, Berg HC (2010) A modular gradient-sensing network for 
chemotaxis in Escherichia coli revealed by responses to time-varying stimuli. 
Mol Syst Biol 6: 382. 

27. Si G, Wu T, Ouyang Q, Tu Y (2012) Pathway-based mean-field model for 
Escherichia coli chemotaxis. Phys Rev Lett 109: 048101. 

28. SpudichJL, Koshland DE, Jr. (1976) Non-gcnctic individuality: chance in the 
single cell. Nature 262: 467-471. 

29. Masson JB, Voismnc G, Wong-Ng J, Celani A, Vergassola M (2012) 
Noninvasive inference of the molecular chemotaetic response using bacterial 
trajectories. Proe Natl Acad Sci U S A 109: 1802-1807. 

30. Berg HC, Tedesco PM (1975) Transient response to chemotaetic stimuli in 
Escherichia coli. Proe Natl Acad Sci U S A 72: 3235-3239. 

31. Park H, Pontius W, Guet CC, Marko JF, Emonet T, et al. (2010) 
Interdependence of behavioural variability and response to small stimuli in 
bacteria. Nature 468: 819-823. 

32. Andrews BW, Yi TM, Iglesias PA (2006) Optimal noise filtering in the 
chemotaetic response of Escherichia coli. PLoS Comput Biol 2: el54. 

33. Clark DA, Grant LC (2005) The bacterial chemotaetic response reflects a 
compromise between transient and steady-state behavior. Proe Natl Acad 
Sci USA 102: 9150-9155. 

34. de Gennes PG (2004) Chemotaxis: the role of internal delays. Eur BiophysJ 33: 
691-693. 

35. Mears PJ, Koirala S, Rao CV, Golding I, Chemla YR (2014) Escherichia coli 
swimming is robust against variations in flagellar number. cLifc 3: e01916. 

36. Schnitzcr MJ (1993) Theory of continuum random walks and application to 
chemotaxis. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 48: 
2553-2568. 

37. Othmcr HG, Xin X, Xuc C (2013) Excitation and Adaptation in Bacteria-a 
Model Signal Transduction System that Controls Taxis and Spatial Pattern 
Formation. International journal of molecular sciences 14: 9205-9248. 

38. Tu Y, Berg HC (2012) Tandem adaptation with a common design in 
Escherichia coli chemotaxis. Journal of molecular biology 423: 782—788. 

39. Zhu X, Si G, Deng N, Ouyang Q, Wu T, et al. (2012) Frequency-Dependent 
Escherichia coli Chemotaxis Behavior. Phys Rev Lett 108: 128101. 

40. Angeli D, Ferrell JE, Jr., Sontag ED (2004) Detection of multistability, 
bifurcations, and hysteresis in a large class of biological positive -feedback 
systems. Proe Natl Acad Sci U S A 101: 1822-1827. 

41. Oleksiuk O, Jakovljevic V, Vladimirov N, Carvalho R, Paster E, et al. (2011) 
Thermal robustness of signaling in bacterial chemotaxis. Cell 145: 312—321. 

42. Turner L, Ryu WS, Berg HC (2000) Real-time imaging of fluorescent flagellar 
filaments. J Bacteriol 182: 2793-2801. 

43. Saragosti J, Calvez V, Bournaveas N, Perthame B, Buguin A, et al. (2011) 
Directional persistence of chemotaetic bacteria in a traveling concentration 
wave. Proe Natl Acad Sci U S A 108: 16235-16240. 

44. Nicolau DV,Jr., Armitagc JP, Maini PK (2009) Directional persistence and the 
optimality of run-and-tumble chemotaxis. Comput Biol Chem 33: 269—274. 

45. LocseiJT (2007) Persistence of direction increases the drift velocity of run and 
tumble chemotaxis. Journal of Mathematical Biology 55: 41—60. 

46. Emonet T, Cluzel P (2008) Relationship between cellular response and 
behavioral variability in bacterial chemotaxis. Proe Natl Acad Sci USA 105: 
3304-3309. 



PLOS Computational Biology | www.ploscompbiol.org 



11 



June 2014 | Volume 10 | Issue 6 | e1003694 



